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Solving realistic quantum systems coupled to an environment is a challenging task. Here we 
develop a hierarchical functional derivative (HFD) approach for efficiently solving the non-Markovian 
quantum trajectories of an open quantum system embedded in a bosonic bath. An explicit expression 
for arbitrary order HFD equation is derived systematically. Moreover, it is found that for an 
analytically solvable model, this hierarchical equation naturally terminates at a given order and 
thus becomes exactly solvable. This HFD approach provides a systematic method to study the 
non-Markovian quantum dynamics of an open system coupled to a bosonic environment. 
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I. INTRODUCTION 

The theory of open quantum systems [1] has received 
great interest because environment-induced effects are 
important in a wide range of research topics such as quan¬ 
tum information [2] and quantum optics [3] . There were 
considerable studies involving environments modeled by 
either bosonic or fermionic baths (see, e.g., Refs. [1, 4- 
13]), as well as structured environments such as spin- 
chain baths [14]. Conventionally, Markov approximation 
has been extensively used because of its simplicity. In¬ 
deed, this is valid only when memory effects of the en¬ 
vironment are negligible. However, this approximation 
becomes invalid when the system-environment coupling 
is strong or when the environment is structured [1, 15]. 
Therefore, non-Markovian environments have to be con¬ 
sidered in explaining new experimental advances in quan¬ 
tum optics. Also, they must be considered in quantum 
information manipulations in which the environmental 
memory is utilized to control the entanglement dynam¬ 
ics [16]. Therefore, it is vital to have a non-Markovian 
description of the system’s dynamics under the influence 
of the memory effects and the backaction of the environ¬ 
ment. Actually, this has long been a challenging task and 
many theoretical approaches have been developed (see, 
e.g., Refs. [1, 4-10, 17-21]). Among these approaches, 
the non-Markovian quantum state diffusion (QSD) [4- 
6] has been proven to be a powerful tool to study the 
quantum dynamics of the system and exact analytical 
results were derived for many interesting systems such as 
dissipative multi-level atoms [8] and quantum Brownian 
motion [4, 22] which was also analytically solved via a 
path-integral approach [17]. Quantum continuous mea¬ 
surement employing the QSD technique was also stud¬ 
ied [23, 24], 

For most realistic open quantum system problems, it 


is almost impossible to find any useful analytical solu¬ 
tions. Therefore, one has to develop numerical methods 
to study the quantum dynamics of the open systems. 
However, the application of the non-Markovian QSD is 
greatly hindered unless one can cast it to a numerically 
implementable time-local form. Recently, two hierarchi¬ 
cal approaches have been proposed; one is the stochastic 
differential equation (SDE) method [9] based on a func¬ 
tional expansion of a system operator, and the other is 
an approach using a hierarchy of pure states (HOPS) [10] 
to calculate a related functional derivative. Apart from 
these two approaches, an earlier attempt at a perturba¬ 
tive solution of the non-Markovian QSD equation was 
based on hierarchical functional derivative (HFD) [7]. 
However, because of its apparent complexity, the hierar¬ 
chical equations presented there were implemented only 
up to the second order where the higher-order terms were 
approximated by a simplified operator. Here we develop 
a systematic and efficient higher-order HFD approach to 
solving the non-Markovian quantum dynamics of an open 
system coupled to a bosonic environment. Remarkably, 
a compact explicit expression for an arbitrary order hi¬ 
erarchical equation is derived. Moreover, it is found that 
for an analytically solvable model, the hierarchical equa¬ 
tion naturally terminates at a given order, so it becomes 
exactly solvable. Thus, our method provides not only an 
approach for efficiently solving the non-Markovian quan¬ 
tum dynamics of an open system, but also a systematic 
method for the exact solution of analytically tractable 
open systems. 

The paper is organized as follows. In Sec. II, we present 
our HFD method for studying the non-Markovian QSD. 
Then, we show the relationship of the HFD method to 
the recently developed SDE and HOPS method in Sec. Ill 
and Sec. IV, respectively. Finally, discussion and conclu¬ 
sion are given in Sec. V. Moreover, in Appendix A, we 
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present a detailed derivation of the HFD equation for 
a general bath correlation function. The mathematical 
relationship between the SDE approach and our HFD 
method is explicitly demostrated in Appendix B. 


II. THE HFD METHOD FOR 

NON-MARKOVIAN QSD 

We study a generic open system with the Hamilto¬ 
nian [4, 5] 

H = H s + ^ [g k Lb\ + g k L^b)}j + ^ Ukb\b k , (1) 

k k 

where H s is the Hamiltonian of the system of interest, L 
is the coupling operator called Lindblad operator, b k is 
the fcth mode annihilation operator of the bosonic bath 
with a frequency uj k and g k denotes the coupling strength 
between the system and the bosonic bath. The bath 
state can be specified by a set of complex numbers {z k } 
labeling the coherent state of all bath modes and the ef¬ 
fect of the bath is characterized by the zero-temperature 
bath correlation function a(t,s) = J2k \g k \ 2 e~ lulk( ' t ~ s \ 
Defining zf = — i J^k 9k z k elulkt > one can interpret z k as 
a Gaussian random variable and z% becomes a Gaussian 
process with its statistical mean given by the bath cor¬ 
relation function a{t, s) = (( ZtZ *)). The wave function 
|= (2*|\E'tot(i)) ! obtained by projecting the quan¬ 
tum state I'FtotW) °f the total system onto the bath state 
\z), is called a quantum trajectory and obeys a linear 
QSD equation 

= [-iH s + Lz*-tfd(t,z)\ I(2) 

where O is an operator defined by the functional deriva- 
tive 31*1 = 0{t,s,z*)\ip z *(t)) and 0(t,z*) = 

fo a(t, s)0(t,s,z*)ds [4]. Non-Markovian master equa¬ 
tions can also be obtained using this approach [6, 22]. 
For this linear QSD, | ip z *(t)) is an unnormalized wave 
function, so it can be of different orders of magnitude at 
different evolution times. This gives rise to an inefficient 
Monte Carlo sampling in numerical calculations. Thus, 
one utilizes a non-linear QSD [4] for the normalized state 

= [-iH„ + A t {L)~z* t - A t {P)0{t, z*) 

(3) 

which is derived via the Girsanov transformation zf = 
z t + fo a *(t> s)(L^) s ds. Here A t (L) = L — (L) t , and the 
reduced density operator p s (t) = Tr en v| v I'tot)('I / tot| can 
be obtained from the ensemble average of the normalized 
quantum trajectories as p s (t) = ((\tp£* (t))(Vtz(Q|)}- 
The non-Markovian QSD equation (3) is exact, but the 
key challenge in solving it relies on the determination of 


the O operator. It is difficult to analytically obtain the 
explicit expression of the O operator except for a few sim¬ 
ple models such as dissipative multi-level atoms [8], the 
quantum Brownian motion [4] and dissipative multiple 
qubits [25]. Owing to the importance of open systems, 
efforts [7, 9, 10] have been devoted to develop numeri¬ 
cal methods for solving Eq. (3). In Ref. [7], an approach 
was proposed to calculate the functional derivative by 
defining a set of Q k operators as 

Q k (t,z*) = J a(t,s)J^Q k -i(t,z*)ds, (4) 

where Qo{t, z*) = 0(t,z*). Because of the complexity 
involved, an explicit hierarchical equation of motion for 
this set of operators was obtained in Ref. [7] only up 
to the second-order Q 2 operator while the higher-order 
terms were approximated by a simplified operator. Below 
we give an explicit hierarchical equation of motion up to 
any high orders and also applicable to an arbitrary bath 
correlation function a{t, s). For this purpose, we include 
the time derivatives of a(t, s ) and define 


Q k \t,z*) = - p(A)'p(A-i) .. yp(A) 


x [ ds 0 a < ' JO \t,s 0 )0(t,so,z*), (5) 


Jo 

where j = J2, Ji e i = (jo, ji, ■ ■ ■,jk ) with e ; being the ith 
unit vector, "pt?) = / Q 4 dsa^\t, s) and = 

■r~a(t,s). It is derived (see Appendix A) that this set of 
operators satisfy 


dt 


Qk\t,z*) = 5>W(0) L, z*) 


+ ^2Q ( i +ei \t,z*)-L^Q^(t,z*) 


i =0 


-iH sys + Lzl Qk(t, z*)\ + 6 0 , k a^(0)L 

EE [tfQf’ Ci \t,z*),Q^4-'\t,z*)\ , (6) 


i=0 Cj 


where Q^\t,z*) = Q k {t,z*), a^(0) = (t,t), 

£>CM) = {jo, ■ ..ji-i,ji+i,... ,j k ) excludes ji in the vec¬ 
tor j, and Ylc- the sum of k\/[i\{k — z)!] terms which 
include all possible cases of choosing i elements from a k- 
dimensional vector {j \,..., j k ). Using functional expan¬ 
sions [6], we can also immediately prove that Q k \t,z*) 

operators have such a symmetry that Q k °’^’"'^ k \t,z*) 
coincide for all permutations of (ji, j 2 , • ■ ■, jk)- This 
enables us to greatly reduce the number of equations 
that should be solved, because we only need to calcu¬ 
late those Q k \t,z*) operators with normal-ordered ji 
(ji < J2 < • • • < jk)- 

It is also worth pointing out that since the QSD equa¬ 
tion for a finite-temperature bath with self-adjoint Lind¬ 
blad operator L (i.e., L = L^) is formally the same as the 
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zero-temperature one with a different finite-temperature 
bath correlation function a(t,s ) [4], so Eq. (6) also ap¬ 
plies to that finite-temperature case. Note that the finite- 
temperature case with a self-adjoint L covers a large num¬ 
ber of open-system problems, including the important 
spin-boson model without the rotating-wave approxima¬ 
tion (see, e.g., [9]). 

When modeling the structure of the bosonic bath, 
an important and widely used choice is the Lorentzian 
spectrum, which corresponds to an environmental noise 
z t of the Ornstein-Uhlenbeck type with autocorrelation 
a(t,s) = T7exp(—y|t — s\)/2. This kind of bosonic 
bath has been used to describe many interesting prob¬ 
lems [1, 3], and it is easy to observe a non-Markovian to 
Markovian crossover by just increasing 7. In addition, 
this can greatly simplify the hierarchical equations since 
oft\t,s) = and Q^ = (—7 ) jB Qk, where 

J s = y~b jj ■ In this case, by utilizing our general re¬ 
sult (6), a compact hierarchical equation for the Q k op¬ 
erators given by Eq. (4) can be explicitly written as 

§- t Qk{t, z*) = MO) [L, Qk-i(t, z*)\ + Va(0)L 
~(k + 1)7 Q k (t, z*) + [-Ws + L%, Q k (t, z *)] (7) 

k 

- L^Q k+1 (t, z*) **), 5*)] , 

i=0 

where Cf = k\/[i\(k — z)!] is the binomial coefficient 
and the initial conditions are Qfc(0,S*) = 0 for k > 
0. This set of hierarchical equations can be numeri¬ 
cally solved perturbatively if terminated at order Af + 1 
by putting either Qjv+i(t,z*) = 0 or Qjy+i{t,z*) = 
fo dsa(t, s) [L, Q/s{t, 5*)] [7]. 

Now we also explore the use of this approach as a sys¬ 
tematic method for models with exact solutions. For 
open systems that are known to be exactly solvable by 
QSD, the number of expansion terms in 

0(t, z*) = O i0 \t) + [ 0 ( - 1 ' > {t,vi)z* 1 dvi 

Jo 

+ [ [ O i2) (t,v 1 ,v 2 )z* i z* 2 dv 1 dv 2 +■■■■ (8) 

J 0 J 0 

must be finite [6]. Thus, O^ = 0 for all n larger than a 
given finite integer N c . In this case, we can readily show 
that 

QN c (t) =N C \ [ ■■■ f ...a(t,v Nc ) 

Jo Jo 

x O iNc \t,v 1,. ,.v Nc )dvi. .. dv Nc , (9) 

which is independent of the noise z $. Therefore, its func¬ 
tional derivative with respect to z* is zero. From Eq. (4), 
it follows that Qk = 0 for all k > N c + 1, so that the hier¬ 
archical equation naturally terminates and the approach 
becomes exact. This provides a useful and systematic 



t 
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FIG. 1. (Color online) (a) Ensemble average values of the an¬ 
gular momentum (Jx), ( J y ), and (J z ) versus time t. (b) En¬ 
semble average values of the trace norm of Qo(t, z*), Qi (t, z*), 
and Qk(t, z*), with k > 2, as a function of t. In both (a) and 
(b), we use 1000 noise realizations and 7 == T = uj = 1. The 
solid curves correspond to the analytic solutions; the solid, 
open circles and triangles correspond to the results obtained 
using our HFD approach. 


method to deal with an open system with unknown prop¬ 
erties by just implementing the hierarchical equation; if 
the hierarchical equation has a natural termination at a 
given order, the considered model is analytically solvable 
and our results will be essentially accurate. 

As an illustrative example, we consider an analytically 
solvable three-level system [8], with H sys = ujJ Zi and 
L = J-. The functional expansion of its O operator 
is only up to the first-order term, i.e., N c = 1. Thus, 
only Qa{t,z*) and Q\[t) are involved in the hierarchical 
equation, and Q k = 0 for k > 2. We numerically solve 
the hierarchical equation (7) up to order J\f = 10. From 
Fig. 1(a), it can be seen that the numerically calculated 
ensemble averages (Jj), z = x, y and z, agree well with 
the exactly solved results. Also, the trace norms 11 Qfc|[ = 

Tr are calculated in Fig. 1(b), which shows 

that for k > 2, the Qk operators remain zero and the 
hierarchical equation naturally terminates at order k = 
2 (i.e., N c = 1), in full consistency with the analytical 
derivations. 
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FIG. 2. (Color online) (a) Ensemble average (a z ) versus time 
t for 7 = 0.2, 0.4 and 0.8, obtained using the SDE approach 
(solid curves) and the HFD method (solid, open circles and 
triangles). Here Af = 30, Ty = 0.2, and 1000 noise realiza¬ 
tions are used, (b) Total number of equations versus termi¬ 
nation order Af, where the curve marked by triangles (solid 
circles) corresponds to the number of coupled equations of mo¬ 
tion that should be simultaneously solved in the SDE (HFD) 
approach. 


proximation (RWA) was studied using this approach by 
considering an Ornstein-Uhlenbeck noise (see [9]). This is 
a typical example where the explicit form of the O oper¬ 
ator is unknown and the hierarchical approach can serve 
as a powerful numerical tool. In Fig. 2 (a), we display 
the expectation value of the angular momentum (a z ) as 
a function of time using both the SDE and our HFD ap¬ 
proaches. Here we also use an Ornstein-Uhlenbeck noise, 
so as to directly compare with the SDE results. An ex¬ 
cellent agreement is reached for the numerical results ob¬ 
tained by them. As a matter of fact, these two hierarchi¬ 
cal methods are very closely related mathematically, and 
we prove in Appendix B that the operators Qk and Qm^ 
are related by 

S4M-) = £ (11) 

n=k ' ' 

The key advantage of the HFD method over the SDE 
approach is that the HFD equation (7) effectively groups 
the operator according to Eq. (11) and sums up their 
contributions. This makes the HFD method much more 
efficient than the SDE approach. In fact, for a given ter¬ 
mination order Af , the SDE approach needs to simultane¬ 
ously solve coupled equations for Qm' 1 where n + m < Af, 
resulting in a total number of i(7V+2)(A/’+l) equations. 
On the other hand, the HFD approach greatly reduces 
the total number of equations to Af + 1. Figure 2(b) 
shows the total number of coupled differential equations 
that should be solved in both the HFD and the SDE ap¬ 
proaches. The very high efficiency of the HFD method 
over the SDE approach is clearly seen when increasing 
Af. 


III. THE RELATIONSHIP TO THE SDE 
METHOD 


IV. THE RELATIONSHIP TO THE HOPS 
METHOD 


Recently, a numerical approach [9] was developed to 
solve the non-Markovian QSD equation via a set of 
stochastic differential equations (SDE). A key part of the 
SDE formulation is the introduction of a Q operator 




vi) v m )z* m+1 


...IJ n O (n) (f,!)i,...,!) n )dj;i ...dv n , ( 10 ) 


where 0^ n \t : vi,...,v n ) corresponds to the_nth-order 
functional expansion term in Eq. (8) and 0(t, z*) = 

Sn=o Qo"^i ■?*). For m ^ 0, Qm' 1 do not directly con¬ 
tribute to 0(t, z*) but form a set of hierarchical equations 
with Qq and need to be solved simultaneously. Within 
this framework, one can calculate the quantum trajectory 
up to an arbitrary order of the environmental noise. 

The spin-boson model without the rotating-wave ap- 


Previously, almost all QSD approaches are focused on 
how to calculate the functional derivative associated with 
the O operator. Recently, a hierarchy of pure states 
(HOPS) approach [10] was developed as a numerical tool 
which, instead of using the O operator, introduces a set 
of pure states 


IV’fcW) = J a(t,s)J^ds\i> k - 1 (t)), (12) 

where \ipo{t)} = |'0(^))- hr this approach, a set of hierar¬ 
chical equations of motion was found for \ipk{t))- Its ad¬ 
vantage is that the hierarchical equations deal with state 
vectors of size diin("H s ) x 1 rather than the operators of 
size dim("H s ) x dim("H s ), where dim("H s ) is the dimension 
of the system’s Hilbert space. From the definition, it is 
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easy to see that 

IVhO)) =Qoit,z*)\ip 0 (t)), 

1^2 (t)} =Qi(t,z*)\ip 0 (t)) + Q 0 (i,5*)|^i(t)}, 

\i>3{t)) =Q 2 {t, z*)\ipo(t)) + 2Qi(t, z*) iV'i(i)} 

+ Qo(M*)hMi)), ., 

and in general, 

fc-i 

i^ = E c E l2 ^ r )i^-- i >- ( 13 ) 

j=0 

Therefore, it is also possible to formulate the HOPS ap¬ 
proach using the Qk operators in Eq. (4). As shown 
above, for an exactly solvable model, Qi(t,z*) = 0 for 
all i > N c + 1, so that the hierarchical equations for 
Qi(t , z*) terminate at the lV c th order. It is interesting to 
note that there are infinite number of hierarchical equa¬ 
tions for \ipk) in the HOPS approach, but only N c + 1 
nonzero operators Qo(t, z*), Qi(t, z*), ..., Q.N c {t, z*) are 
needed to express all of these \ipk)- This reflects the exact 
solvability of the considered model. 


V. DISCUSSION AND CONCLUSION 


naturally terminates at a given order and becomes ex¬ 
actly solvable for an analytically solvable model, it pro¬ 
vides a systematic perturbation for a generic open system 
irrespective of the existence of the time-local O operator. 
Our HFD method applies to an arbitrary bath correlation 
function. Also, it can be naturally extended to the case of 
a finite-temperature bath when the Lindblad operator in 
the interaction Hamiltonian is self-adjoint, including the 
important spin-boson model without the rotating-wave 
approximation. 
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Appendix A: Derivation for the HFD equation with 
a general bath correlation function 


Solving non-Markovian dynamics of an open quantum 
system has long been a challenge. The conventional 
master equation for the system’s reduced density ma¬ 
trix is driven by a super-operator 1C, which is given by 
the perturbation expansion with respect to the strength 
of system-bath interaction [1], Whereas the expansion 
could be done manually up to the orders beyond 2, it 
is difficult to find an automatic numerical way to obtain 
the higher-order 1C operator. As such, the conventional 
master equation is used in the weakly-coupled scenario. 
Likewise, the quantum state diffusion equation driven by 
an O operator (which plays a role similar to 1C) encoun¬ 
ters the same problem. As a breakthrough, this work 
develops a systematic and efficient higher-order HFD ap¬ 
proach for solving the non-Markovian quantum trajecto¬ 
ries of an open system coupled to a bosonic environment. 
A compact explicit expression for an arbitrary order hi¬ 
erarchical equation of motion is derived and it can be 
efficiently implemented numerically. As a distinctive ad¬ 
vantage of this method, while this hierarchical equation 

I 


In this appendix, we derive the HFD equation [i.e., 
Eq. (6)] for an open system with a general bath correla¬ 
tion function a(t, s). Define 

Qk\t,z*) =pCik)'pUk-i) -pCh) 

x f ds 0 a ( ' J °\t,s 0 )O(t,s 0 ,z*), (Al) 
do 

where j = Y,,ji e i = (jo, ji, • • ■ > jk) with e ; being the ith 
unit vector, = f* dsa^ (t, s)j^r, and a®(f,s) = 

ja(t,s ). At the zeroth order, using the quantum state 
diffusion (QSD) equation 


0(t,s,z*)] = 


—iH sys + Lz* - z*), 0(t, s, z*) 

(A2) 


rf S(4 0) (t,z*) 

5z* 


we have 


d n Uo)r - 9 


dt 


(*>**) = 


dt 


ds 0 a M (t, s 0 )O(t, s 0 , z*) 


= a (jo) (0 )0(t,t,z*) + J ds 0 -^ (a (jo) (t, s 0 )) 0(f, s 0 , z*) + J ds 0 a M (t,s)6(t,s 0 , z*) 

= a^{0)L + Q <9o+1 \t,z*)-L^Q[ 0do \t,z*)+\-iH sys + Lz:-L^Q { °\t,z*),Q { 0 jo \t,z*)} ■ (A3) 
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For fc > 1, it can be seen that 


IpU) = 1 
dt dt 


UO 


= a W) (0) 


dsa^\t, s) 
<5 


5z* 


5zl 

dsa^ +1 \t, s ) 


S J 

* 5 


Szt 


= a {, ' ) (0)A+P Ci+1) - 


5z* 


(A4) 


Here we define 2?(j ,i) = {jo, ... j*_i, ji+i,..., jfe), which 
excludes ji in the vector j. For any n-dimensional 
vector v = (vi,...,v n ) and integer x, the notation 
(x, v) represents a new (n + l)-dimensional vector v' = 
(x,vi,..., v n ). We also introduce a vector Ci that chooses 

I 


i elements from {ji,... ,jk) and its complement Cj that 
contains the remaining elements. For example, for 
(1,2,3,4, 5), C2 can be (1,4) and then C2 = (2, 3,5). 

We thus have 




dt 

k Q'pUi) 


— ptik)ptik- 1) . 


dt 


. pdd g0’°) 


. 'p(ik)'pUk- 1) -p(ii) d2o°V>z*) 

at 


(A5) 


Using Eq. (A4) for the first term on the RHS of Eq. (A5) 
and Eq. (A3) for the second term, we have 


dt 


Q,f{t,z*) =^>(A)p(A-i) ...p(h+i) a a«)(o)Ap(ii-i). ..pWQU°\t, Z *) 

0 Zj. 


2=1 


Note that 


+ 'pUk)'p(jk-i)' < /pOi+i) _ _ < 'p( <?1 )Qo 7 °)(t, z*) 


2=1 








fc+i v c 5 z / 


_ V Uk) v Uk- 1).. ,p(A) [Lt Q (0)( tjZ *) ) qCj°) (*,*•) 


(A6) 


pdk) v Uk-d _ pUi+ 1 ) _ v (h) qUo)^ z *^ = g0+ e ‘)(t, ar*), 

pCA)p(A-i) _ _. p(h+i) Q Oi)(o)^p(^-i)... V {jl) Q ( 0 jo) (t, z*) = a Ui \0) 

OZj 


L,o£H A) (t,z*) 


Then, we finally have 


(A7) 


= E tt(ji) (°) [ L , + E 2i j+ei) (M*) + Qi j+eo) (M*) - A t Qfc°| J i(t,^*) 


2=1 


+ [-zi7 sys + L<,Q^ j) (t,z*)] -vtik) v Uk-i)'" V Ui) ^Q^\t,z*),Q^°\t,z*) 

E« (ii) (0) +E2i j+ei) (i^*)-^ t Qi+i ) (i^*) 

2=1 2=0 

zij sys + L4,g( j) (t,z*)] - EE 


2=0 Ci 


(AS) 


where ]T) c is the sum of C k = k\/[i\{k — *)!] terms that 
include all possible cases of choosing i elements from a 
fc-dimensional vector {j\,, jk)- This comes from the 
term 


v Uk) v Uk-i) ...-p(A) \^Q^{t,z*),Q^ 0 \t,z*)\ , (A9) 


where we have chosen i operators in pdk)p(jk-i) _ _ pU i) 
to act on QQ°\t,z*), and the remaining (fc — i) 
•p(A) operators act on Qo°\t, z*). For example, 
for pUs)pUi)pUa)pU3)pUi) } one possible case with 
i = 2 is that pUa)pih) ac ts on Qo°\t,z*) to give 
Q ( 2 °’ jlj3) (t, z*), and pUdp(n)p(h) acts on Q^°\t, z *) to 



























7 


give Q^ 0 ’- 72 > J4J5 )(i ) z *)' For the operator, termina¬ 
tion of the index k is determined by the functional deriva¬ 
tive of the O operator, whereas termination of the index 
j is determined by the bath correlation function a(t,s). 
In numerical calculations, this set of equations can be 
terminated by imposing Q^ 1 = 0 for all k + Yi ji > A/", 


where A/" is a suitably chosen integer. 

As an example, we use an exactly solvable model to 
show what the HFD equations look like in the pres¬ 
ence of a general bath correlation function. Consider 
a three-level atom in a multi-mode bosonic bath [8], 
where H sys = ojJ z , and L = J^. The HFD equation 

for Q^ j) (i, z*), with j = (j 0 ), is given by 


§- t Qo o) (t > z *) = a (jo) (0 )L + Qo° +1 \t, z*) - L^Q^ j0 \t,z*) + [-*ff sys + Lz* - z*), Q ( 0 30 \t, z*) 


, (A10) 


where jo = 0,1,2,.... The HFD equation for Q^\t, z*), with j = (jo, ji), is given by 


d 


f t Q[ 30 ’ n) (t,z*) =«<*>(0) [L,Qtf°\t,z*) 


Q[ j0+1 ’ jl \t, z*) + Q[ J0 ’ n+L, {t, z*) + [~iH sys + Lz*, 
^Q^\t,z*),Q[ JO ’ n \t,z*)] - \^Q^ jo \t,z*),Q ( 0 j °\t,z*) 


l(j'o iJi + 1)/ 


UoJi)^z*) 

(All) 


where jo, ji = 0,1, 2,... and jo < ji- Because the func¬ 
tional expansion of the O operator terminates at the first 
order for this exactly solvable three-level model, it follows 
from Eqs. (Al) and ( 8 ) that 

Q<j\t,z*) = 0, j = (j 0 , ji, - - -, j fc ), k> 2 . (A 12 ) 

Solving HFD equations in Eqs. (A10) and (All), we 
can obtain O for an arbitrary bath correlation function 
a(t, s ) and then determine the quantum-dynamical be¬ 
havior of the system using the QSD equation. To ter¬ 
minate Eqs. (A10) and (All) in numerical calculations, 
we can impose g[ l im< “ +1) = a Qjf max) and Qp= 
a Qi m “ Jm “), where j max is a suitably chosen integer and 
a is a parameter determined by the bath correlation func¬ 
tion. 

In particular, for an Ornstein-Uhlenbeck noise, we have 

Qo 1 ' = ~72 o 0) = ^7So and qJ 0,1) = - 7 Qi°’ 0) = - 7 Q 1 . 
Thus, j max = 0 and a = — 7 . Then, Eqs. (A10) and (All) 
simply reduce to 

lQ 0 (t,z*)=a(0)L-'yQ 0 (t,z*)-L'Q 1 (t,z*) 

+ [-iH sys + Lzt-L'Q 0 (t,z*),Q 0 (t,z*)\ , (A13) 

and 

§- t Qi(t, z*) = a(0) [L, Q 0 (i, z*)\ - 2 7 Q 1 (t, z*) 

+ [~iH sys + Lz*, Qi{t,z*)] - [L'Q 0 (t,z*), Qi{t,z*)\ 
-[L'Qi&z^Qo&z'j], (A14) 


Appendix B: The relationship between SDE and 

HFD methods 

In this appendix, we explicitly show the relationship 
between the stochastic differential equation (SDE) ap¬ 
proach [9] and our generalized hierarchical functional 
derivative (HFD) approach. 

From the definition of the Qm 1 operator in Eq. (10), 
we have 

J q dsa{t,s)-^Q ( £\t,z*) = (n - m)Q ( n n J +1 . (Bl) 
Thus, from Eqs. (Al) and (4), it follows that 

Q 0 (t,z*) = d(t,z*) = Z2Q£ l) (t,2*), 

n =0 

Qi(t,z*)= J dsa(t,s)J^Q 0 (t,z*) 

JV 

= 5>Q[ n) , (B2) 

n— 1 

Q2(t,z*)=J dsa(t,s)-^Q 1 (t,z*) 

= ^2n{n-l)Q^\ 

n =2 


By mathematical induction, it can be obtained that if 
a r 

Qk(t, z*) = ^2 n ( n - 1) • ■ • (n - k + 1 )Q ( j! 1 ' > 

n—k 

= (M-), <B3) 


which give the results in Ref. [8]. 
















then for Q k +i, we have 


This proves the result given in Eq. (11). 


Q k+1 {t, z*) = J2 (n J q dsa(t, 5*) 


n—k 

a r 


= E 


(n — k — (®4) hierarchical equation of motion for the Q^' 1 operator is 

n=k+i derived in Ref. [9] as 

I 


When an Ornstein-Uhlcnbeck noise is considered, the 


d ^(n) / ■ ~*\ ~ ^ 


^ 7 Qk (t, z*) =5 nfi a(0)L + 


dt 


max{l, n} 


«(0) \L,Qt~i\t,z*) 


(n-k) 
max{l, n} 




-(fc + l) 7 Q< n) (i,r) + l-iH^Q^^z*)] -(n+l^Q^^z*) 


-EE 


s~ip s~in—p 
° n-k-l 


fin 

p—0 l k 




p—l\ u ->* J’> y °6k—p+l\ 
_ sr^Af n \ d r\i n ) ( 


(B5) 


Substituting Eq. (B5) into ^ Q k (t,z *) = J2n=k (n-k) 1 . (t,-5*), h f°U° ws that for k = 0, one has 

a r 


9 Qo{t,z*) =J2 ^:Qo l) {t, z*^ 


dt 


n —0 


a r 


n =0 


=a(0)L + J2{% L, Qo ll \t, z*) - 7 Q<j n) (t, 5*) + -tH sysl Q { 0 n \t, z*) 


n p ftp fin—p 

1 n ~ l -(n + l)lSQ<T +1) (t,z*) 


EE-- ^'^p-lVi* h^i-p Vi*) -yf-r vr 

p —0 l=p o 

=a(0)L + ~z* t [L, Qo(t, 2*)] - 7 Q 0 (i, 5*) + [~iH sys , Q„(*, 5*)] - [L*Q 0 (t, z*), Qo(t, 5*)] - L'Q^t, ~z*), 

(B 6 ) 


and for k > 1 , one has 


d 


a r 


dt Qk A 5 *) -E (n _ fc)! a ^fc n E 5 *) 

n! f k 
(n — k)\ \ n 


%=k 

a r 


= £t^M *“«»[■ 


2=k 


(n — k) 


l,q { k 1) (t, ^*)1 - (fc + i) 7 ^ n) (t,r) 


+ 


n y^pfiin—p 

—iH sys , z*)] - (n + l)I' + < 5 fc"E(t, **) - E E 7 jn' t ~ i 5 *), r 


p=0 l k 

=ka( 0) [I/, Q fc _i(t, «*)] + z* t [L, Qk(t, 5*)] — (k+ 1 ) 7 Q fc (i, 5*) - L+Qk+i^, 5*) 


AT 


[ *i?sys, Qfc(t, 5*)] X] , _ ,w EE 


s~ip s~in—p 
°Z ^n-k-l 


(n — k)\ 

n=k v ' p —0 / 

The last term on the RHS of Eq. (B7) can be rewritten as 


CV 


L'Q^[tX),Qt7U^ 


(B7) 


A/ - 


y —-—y 

*-^ (n — k)! z— 1 


n—k 


fPfin-p 


p,Z 


CT? 


= EE 


k\ 


(p - l)!(k + 1 —p)\ 


n—k p,l 

= Y. C i [L'Qi(t,z*),Qk-i(t,z*)\ , 


it^lnCp) (t z*) ( n p ) ! n (ra ~ p ) (t z*) 
l\^P-^ Z >' ( n -k-l)\ Qk -v +A ' j 


(B8) 
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where the substitution p — l —> i is taken in the last line. 
Finally, substituting Eq. (B 8 ) into Eq. (B7), we obtain 


Eq. (7), i.e., the hierarchical equation of motion for the 
Qk operator. This proves that Eq. (B5) is mathemati¬ 
cally equivalent to Eq. (7). 
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